home *** CD-ROM | disk | FTP | other *** search
/ The Atari Compendium / The Atari Compendium (Toad Computers) (1994).iso / files / prgtools / gnustuff / tos / updates / update15.zoo / pml / pmlsrc / log.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-02-04  |  9.9 KB  |  431 lines
  1. /************************************************************************
  2.  *                                                                      *
  3.  * A replacement log() routine for PML library.  It really computes     *
  4.  * base 2 logarithm which can be multiplied by a proper constant        *
  5.  * to provide final answer.  A rational approximation of an original    *
  6.  * is replaced by a polynomial one.  In practice, with a software       *
  7.  * floating point, this gives a better precision.                       *
  8.  *                                                                      *
  9.  * Michal Jaegermann, December 1990                                     *
  10.  *                                                                      *
  11.  * If __GCC_HACK__ is defined then we are folding log and log10 routines*
  12.  * by making in assembler an extra entry point.  Do not define that     *
  13.  * for portable routine!!                                               *
  14.  *                                                                      *
  15.  * 68881 support added by Michael Ritzert, November 1991        *
  16.  ************************************************************************
  17.  */
  18.  
  19. #ifdef __GNUC__
  20. #  ifndef __GCC_HACK__
  21. #    define __GCC_HACK__
  22. #  endif
  23. #  if !(defined(atarist) || defined(atariminix) || defined(mc68000))
  24. #    undef __GCC_HACK__
  25. #  endif
  26. #endif
  27.  
  28. /************************************************************************
  29.  *                                    *
  30.  *                N O T I C E                *
  31.  *                                    *
  32.  *            Copyright Abandoned, 1987, Fred Fish        *
  33.  *                                    *
  34.  *    This previously copyrighted work has been placed into the    *
  35.  *    public domain by the author (Fred Fish) and may be freely used    *
  36.  *    for any purpose, private or commercial.  I would appreciate    *
  37.  *    it, as a courtesy, if this notice is left in all copies and    *
  38.  *    derivative works.  Thank you, and enjoy...            *
  39.  *                                    *
  40.  *    The author makes no warranty of any kind with respect to this    *
  41.  *    product and explicitly disclaims any implied warranties of    *
  42.  *    merchantability or fitness for any particular purpose.        *
  43.  *                                    *
  44.  ************************************************************************
  45.  */
  46.  
  47.  
  48. /*
  49.  *  FUNCTION
  50.  *
  51.  *    log   double precision natural log
  52.  *
  53.  *  KEY WORDS
  54.  *
  55.  *    log
  56.  *    machine independent routines
  57.  *    math libraries
  58.  *
  59.  *  DESCRIPTION
  60.  *
  61.  *    Returns double precision natural log of double precision
  62.  *    floating point argument.
  63.  *
  64.  *  USAGE
  65.  *
  66.  *    double log (x)
  67.  *    double x;
  68.  *
  69.  *  REFERENCES
  70.  *
  71.  *    Computer Approximations, J.F. Hart et al, John Wiley & Sons,
  72.  *    1968, pp. 105-111.
  73.  *
  74.  *  RESTRICTIONS
  75.  *
  76.  *    The absolute error in the approximating rational function is
  77.  *    10**(-19.38).  Note that contrary to DEC's assertion
  78.  *    in the F4P user's guide, the error is absolute and not
  79.  *    relative.
  80.  *      ** modified ** theoretical precision is only 10**(-16.5);
  81.  *      it works better in practice.
  82.  *
  83.  *    This error bound assumes exact arithmetic
  84.  *    in the polynomial evaluation.  Additional rounding and
  85.  *    truncation errors may occur as the argument is reduced
  86.  *    to the range over which the polynomial approximation
  87.  *    is valid, and as the polynomial is evaluated using
  88.  *    finite-precision arithmetic.
  89.  *    
  90.  *  PROGRAMMER
  91.  *
  92.  *    Fred Fish
  93.  *      Modifications - Michal Jaegermann
  94.  *
  95.  *  INTERNALS
  96.  *
  97.  *    Computes log(X) from:
  98.  *
  99.  *      (1)    If argument is zero then flag an error
  100.  *        and return minus infinity (or rather its
  101.  *        machine representation).
  102.  *
  103.  *      (2)    If argument is negative then flag an
  104.  *        error and return minus infinity.
  105.  *
  106.  *      (3)    Given that x = m * 2**k then extract
  107.  *        mantissa m and exponent k.
  108.  *        (3a)  If m = 0.5 then we know what is the final
  109.  *              result and calculations can be shortened.
  110.  *
  111.  *      (4)    Transform m which is in range [0.5,1.0]
  112.  *        to range [1/sqrt(2), sqrt(2)] by:
  113.  *
  114.  *            s = m * sqrt(2)
  115.  *
  116.  *      (5)    Compute z = (s - 1) / (s + 1)
  117.  *        (5a)  Modified - combine steps (4) and (5) computing
  118.  *              z = (m - 1/sqrt(2)) / (m + 1/sqrt(2))
  119.  *
  120.  *      (6)    Now use the approximation from HART
  121.  *        page 111 to find log(s):
  122.  *
  123.  *        log(s) = z * ( P(z**2) / Q(z**2) )
  124.  *
  125.  *        Where:
  126.  *
  127.  *        P(z**2) =  SUM [ Pj * (z**(2*j)) ]
  128.  *        over j = {0,1,2,3}
  129.  *
  130.  *        Q(z**2) =  SUM [ Qj * (z**(2*j)) ]
  131.  *        over j = {0,1,2,3}
  132.  *
  133.  *        P0 =  -0.240139179559210509868484e2
  134.  *        P1 =  0.30957292821537650062264e2
  135.  *        P2 =  -0.96376909336868659324e1
  136.  *        P3 =  0.4210873712179797145
  137.  *        Q0 =  -0.120069589779605254717525e2
  138.  *        Q1 =  0.19480966070088973051623e2
  139.  *        Q2 =  -0.89111090279378312337e1
  140.  *        Q3 =  1.0000
  141.  *
  142.  *        (coefficients from HART table #2705 pg 229)
  143.  *      (6a)    Modified - compute, as a polynomial of z, an
  144.  *              approximation of log2(m * sqrt(2)). Coefficients
  145.  *              for this polynomial are not given explicitely in HART
  146.  *              but can be computed from table #2665, for example.
  147.  *
  148.  *      (7)    Finally, compute log(x) from:
  149.  *
  150.  *        log(x) = (k * log(2)) - log(sqrt(2)) + log(s)
  151.  *
  152.  *        (7a)  log2(x) = k - 1/2 + log(m * sqrt(2)).  Multiply
  153.  *              by a constant to adjust a logarithm base.
  154.  *
  155.  */
  156.  
  157. #if !defined (__M68881__) && !defined (sfp004)
  158.  
  159. #include <stdio.h>
  160. #include <math.h>
  161. #include "pml.h"
  162.  
  163.             /* mjr++                */
  164.             /* use a different routine instead    */
  165.             /* see end of listing            */
  166.  
  167. # define HALFSQRT2        0.70710678118654752440
  168.  
  169. static double log_coeffs[] = {
  170.     2.88539008177793058026,
  171.     0.9617966939212138784,
  172.     0.577078018095541030,
  173.     0.4121983030496934,
  174.     0.32062199711811,
  175.     0.2612917312170,
  176.     0.24451102108
  177. };
  178.  
  179. #ifdef __GCC_HACK__
  180. static double log_of2[] = {
  181.     LN2, 0.30102999566398119521 /* log10(2) */
  182. };
  183. #endif
  184.  
  185. static char funcname[] = "log";
  186.  
  187. #ifdef __GCC_HACK__
  188. /*
  189.  * This fake function header may change even from version to a version
  190.  * of gcc compiler.  Ensure that first four assembler instructions in
  191.  * log and log10 are the same!
  192.  */
  193. __asm__("
  194.     .text
  195.     .even
  196.     .globl _log10
  197. _log10:");
  198. #ifdef __MSHORT__
  199. __asm__("
  200.     link a6,#-24
  201.     moveml #0x3e30,sp@-");
  202. #else
  203. __asm__("
  204.     link a6,#-36
  205.     moveml #0x3f30,sp@-");
  206. #endif  /* __MSHORT__ */
  207. __asm__("
  208.     movel a6@(8),d2
  209.     movel a6@(12),d3
  210.     moveq #1,d6
  211.     jra   lgentry");
  212. #endif  /* __GCC_HACK__ */
  213.  
  214. double log (x)
  215. double x;
  216. {
  217.     int k;
  218.     double s;
  219.     struct exception xcpt;
  220.     char * xcptstr;
  221. #ifdef __GCC_HACK__
  222.     int index;
  223.  
  224.     index = 0;
  225. __asm__("
  226. lgentry:");
  227. #endif
  228.  
  229.     if (x <= 0.0) {
  230.     xcpt.name = funcname;
  231.     xcpt.arg1 = x;
  232.         if (x == 0.0) {
  233.         xcpt.retval = -MAXDOUBLE;
  234.         xcpt.type = SING;
  235.         xcptstr = "SINGULARITY";
  236.     }
  237.     else {
  238.         xcpt.retval = MAXDOUBLE;
  239.         xcpt.type = DOMAIN;
  240.         xcptstr = "DOMAIN";
  241.     }
  242.     if (!matherr (&xcpt)) {
  243.         fprintf (stderr, "%s: %s error\n", xcptstr, funcname);
  244.         errno = EDOM;
  245.     }
  246.     }
  247.     else {
  248.     s = frexp (x, &k);
  249.     if (0.5 == s ) {
  250.         s = -1.0;
  251.     }
  252.     else {
  253.         /* range reduction with s multiplied by SQRT2 */
  254.         s = (s - HALFSQRT2) / (s + HALFSQRT2);
  255.         /* polynomial approximation */
  256.         s *= poly ((sizeof(log_coeffs)/sizeof(double)) - 1,
  257.                    log_coeffs, s * s);
  258.         /* and subtract log2 of SQRT2 */
  259.         s -= 0.5;
  260.     }
  261.     /* add the original binary exponent */
  262.     s += k;
  263.     /* multiply to get a requested logarithm */
  264. #ifdef __GCC_HACK__
  265.     xcpt.retval = s * log_of2[index];
  266. #else
  267.     xcpt.retval = s * LN2;
  268. #endif
  269.     }
  270.     return (xcpt.retval);
  271. }
  272. #endif /* !__M68881__ && ! sfp004    */
  273.  
  274. /*****************************************************************/
  275.  
  276. #ifdef    sfp004
  277.  
  278. __asm("
  279.  
  280. comm =     -6
  281. resp =    -16
  282. zahl =      0
  283.  
  284. ");    /* end asm    */
  285.  
  286. #endif    sfp004
  287. #if defined (__M68881__) || defined (sfp004)
  288.     __asm(".text; .even");
  289.  
  290. # ifdef    ERROR_CHECK
  291.  
  292.     __asm("
  293.  
  294. _Overflow:
  295.     .ascii \"OVERFLOW\\0\"
  296. _Domain:
  297.     .ascii \"DOMAIN\\0\"
  298. _Error_String:
  299.     .ascii \"log: %s error\\n\\0\"
  300. .even
  301.  
  302. | pml compatible log
  303. | m.ritzert 7.12.1991
  304. | ritzert@dfg.dbp.de
  305. |
  306. |    /* NAN  = {7fffffff,ffffffff}        */
  307. |    /* +Inf = {7ff00000,00000000}        */
  308. |    /* -Inf = {fff00000,00000000}        */
  309. |    /* MAX_D= {7fee42d1,30773b76}        */
  310. |    /* MIN_D= {ffee42d1,30773b76}        */
  311.  
  312. .even
  313. double_max:
  314.     .long    0x7fee42d1
  315.     .long    0x30273b76
  316. double_min:
  317.     .long    0xffee42d1
  318.     .long    0x30273b76
  319. NaN:
  320.     .long    0x7fffffff
  321.     .long    0xffffffff
  322. p_Inf:
  323.     .long    0x7ff00000
  324.     .long    0x00000000
  325. m_Inf:
  326.     .long    0xfff00000
  327.     .long    0x00000000
  328. ");
  329. # endif    ERROR_CHECK
  330.  
  331.     __asm("
  332. .even
  333.     .globl _log
  334. _log:
  335.     ");    /* end asm    */
  336.  
  337. #endif    /* __M68881__ || sfp004    */
  338. #ifdef    __M68881__
  339.  
  340.     __asm("
  341.     flognd    a7@(4), fp0    | log
  342.     fmoved    fp0,a7@-    | push result
  343.     moveml    a7@+,d0-d1    | return_value
  344.     ");    /* end asm    */
  345.  
  346. #endif    __M68881__
  347. #ifdef    sfp004
  348.     __asm("
  349.     lea    0xfffa50,a0
  350.     movew    #0x5414,a0@(comm)    | specify function
  351.     cmpiw    #0x8900,a0@(resp)    | check
  352.     movel    a7@(4),a0@        | load arg_hi
  353.     movel    a7@(8),a0@        | load arg_low
  354.     movew    #0x7400,a0@(comm)    | result to d0
  355.     .long    0x0c688900, 0xfff067f8    | wait
  356.     movel    a0@,d0
  357.     movel    a0@,d1
  358.     ");    /* end asm    */
  359.  
  360. #endif    sfp004
  361. #if defined (__M68881__) || defined (sfp004)
  362. # ifdef    ERROR_CHECK
  363.     __asm("
  364.     lea    double_max,a0    |
  365.     swap    d0        | exponent into lower word
  366.     cmpw    a0@(16),d0    | == NaN ?
  367.     beq    error_nan    |
  368.     cmpw    a0@(24),d0    | == + Infinity ?
  369.     beq    error_plus    |
  370.     cmpw    a0@(32),d0    | == - Infinity ?
  371.     beq    error_minus    |
  372.     swap    d0        | result ok,
  373.     rts            | restore d0
  374. ");
  375. #ifndef    __MSHORT__
  376. __asm("
  377. error_minus:
  378.     swap    d0
  379.     moveml    d0-d1,a7@-
  380.     movel    #63,_errno    | errno = ERANGE
  381.     pea    _Overflow    | for printf
  382.     bra    error_exit    |
  383. error_plus:
  384.     swap    d0
  385.     moveml    d0-d1,a7@-
  386.     movel    #63,_errno    | NAN => errno = EDOM
  387.     pea    _Overflow    | for printf
  388.     bra    error_exit    |
  389. error_nan:
  390.     moveml    a0@(24),d0-d1    | result = +inf
  391.     moveml    d0-d1,a7@-
  392.     movel    #62,_errno    | NAN => errno = EDOM
  393.     pea    _Domain        | for printf
  394. ");
  395. #else    __MSHORT__
  396. __asm("
  397. error_minus:
  398.     swap    d0
  399.     moveml    d0-d1,a7@-
  400.     movew    #63,_errno    | errno = ERANGE
  401.     pea    _Overflow    | for printf
  402.     bra    error_exit    |
  403. error_plus:
  404.     swap    d0
  405.     moveml    d0-d1,a7@-
  406.     movew    #63,_errno    | NAN => errno = EDOM
  407.     pea    _Overflow    | for printf
  408.     bra    error_exit    |
  409. error_nan:
  410.     moveml    a0@(24),d0-d1    | result = +inf
  411.     moveml    d0-d1,a7@-
  412.     movew    #62,_errno    | NAN => errno = EDOM
  413.     pea    _Domain        | for printf
  414. ");
  415. #endif    __MSHORT__
  416. __asm("
  417. error_exit:
  418.     pea    _Error_String    |
  419.     pea    __iob+52    |
  420.     jbsr    _fprintf    |
  421.     addl    #12,a7        |
  422.     moveml    a7@+,d0-d1
  423.     rts
  424.     ");
  425. # else    ERROR_CHECK
  426.  
  427. __asm("rts");
  428.  
  429. # endif    ERROR_CHECK
  430. #endif /* __M68881__ || sfp004    */
  431.